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der by using body rotation is revisited in an adjoint-based optimal control frame¬ 
work. The cylinder’s unsteady and fully unconstrained rotation rate is optimized at 
Reynolds numbers between 75 and 200 and over horizons that are longer than in 
previous studies, where they are typically of the order of a vortex shedding period or 
shorter. In the best conhguration, the drag is reduced by 19%, the vortex shedding 
is effectively suppressed, and this low drag state is maintained with minimal cylinder 
rotation after transients. Unlike open-loop control, the optimal control is shown to 
maintain a specihc phase relationship between the actuation and the shedding in 
order to stabilize the wake. A comparison is also given between the performance of 
optimizations for different Reynolds numbers, cost functions, and horizon lengths. It 
is shown that the long horizons used are necessary in order to stabilize the vortex 
shedding efficiently. 
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I. INTRODUCTION 


The flow past a circular cylinder is often used to test and compare the ability of control 
methods to suppress flow instabilities. When the Reynolds number of the flow {Re = 
UooD/u, where Uoo is the inflow velocity, D is the cylinder diameter and z/ is the kinematic 
viscosity) reaches about 47, the interaction between the two shear layers becomes linearly 
unstable, resulting in a supercritical Hopf bifurcation^. For higher Reynolds numbers, a 
stable limit cycle occurs in the form of the well-known Von Karman vortex street, where 
vortices are shed in turn from each side of the body. 

The suppression of the wake instability for Re > 47 has been thoroughly studied both for 
the simplicity of the cylinder geometry and for the industrial interest in reducing bluff-body 
form drag and wake fluctuations. A plethora of control methods have been applied to this 
setup and some of the most important and successful approaches were reviewed by Choi^. 
Usually, flow control strategies are categorized into passive methods, which do not require 
an energy input - such as splitter plates^ or control cylinders^i^ - and active methods, which 
do - e.g. cross-flow body displacement^ or cylinder rotation?^. 

More recently, there has been an increasing interest in closed-loop active control methods, 
whereby the actuation is modified in real-time, based on information from the flow-held. This 
allows set-point tracking, as well as disturbance rejection, and crucially, it has the potential 
to stabilize a system about an unstable operating point, thus only requiring a small amount 
of energy input after transients as shown in several studies^”—. 

The feedback control approach chosen here is adjoint-based optimal control. This method 
consists in finding and applying the control input that locally minimizes a cost function over 
a given time horizon. Although this method can theoretically be implemented as a real-time 
feedback control strategy for some systems, this is not currently computationally feasible for 
fluid flows. For such systems, it can therefore be seen as an optimization technique where 
the control is optimized offline. A key advantage of this method is that it does not require 
an a priori knowledge of the physical mechanism that will minimize the cost function, so 
it can be used to gain intuition about unsuspected but effective control strategies. For 
instance Joe et a/.— applied optimal control to the flow over a separated flat plate with a 
jet located near the trailing edge and found that a phase-locked, pulse-like waveform was 
optimal as it interacted with the vortices shed in the wake in a more efficient and robust 
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way than sinusoidal forcing. Based on these results, Joe et ai— then designed a simple 
phase-locked controller using a similar waveform, leading to a near-optimal performance, 
without the need for expensive estimators. Additionally, the cost of hnding this optimal 
control does not increase regardless of the number of control inputs. For instance Bewley 
et al.— applied optimal control to re-laminarize a turbulent channel flow, where every point 
in the channel was used for blowing and suction. Finally, adjoint-based optimal control can 
be applied to steady and unsteady flows, as well as linear and nonlinear flows. Most other 
control approaches are not this versatile especially with regard to the nonlinearity of the 
Navier-Stokes equations. In many contexts, this approach can therefore be seen as providing 
an upper limit on the achievable performance for a given control conhguration. However, 
because in most cases the problem is nonlinear and non-convex, many studies have found 
that the results can be strongly influenced by the dehnition of the cost functionl^J^ and 
by the length of the optimization horizonl^J^J^. The importance of carefully choosing the 
control problem dehnition is therefore also investigated in this article. 

In the present study, rotation of the cylinder about its main axis is chosen as the control 
method. The experimental work of Tokumaru and Dimotakis^ demonstrated the potential 
for this control method to affect the how-held and reduce the drag of a cylinder by up to 
roughly 80% at Re = 1.5 x 10^ and has prompted many other authors to investigate the ehect 
of oscillatory cylinder rotation on the development of the vortex street at diherent Reynolds 
numbers, e.g.— . In these studies, the cylinder’s sinusoidal rotation needs to be vigorous 
for the control to successfully reduce drag: typically the amplitude of the tangential velocity 
on the cylinder surface is several times larger than the inhow velocity, and the Strouhal 
number associated with the cylinder rotation is of the order of St = fD/Uoo = 1, where 
/ is the dimensional frequency. The open-loop rotation of the cylinder needs to be this 
powerful in order to change the wake’s stability properties: the train of small vortices that 
are generated in each shear layer by the rotation and advected by the free-stream enhances 
the wake symmetry and suppresses the interaction between the shear layers. The associated 
form drag is therefore reduced signihcantly, but unsurprisingly, this actuation method uses 
more energy than it savesl^i^i. 

He et al.— were the hrst to consider the optimal control of a cylinder wake at low Reynolds 
numbers using rotary oscillation. In their study, the optimal control was evaluated starting 
from the best open-loop sinusoidal conhguration {(f) = 3, St = 0.74 for Re = 200, where 
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0 is the amplitude of the sinusoidally oscillating tangential surface velocity). They chose 
to hnd the periodic control waveform (with a several frequency components) that would 
minimize drag, based on a horizon of 3 (high-frequency) forcing periods, and a Reynolds 
number of 200 and 1000. The resulting waveform is qualitatively similar to the open-loop 
one (high amplitude and high frequency), but with marginally improved drag reduction. 
Homescu et alM used a similar setup to minimize the difference between the flow state 
and a target flow-held, which was chosen to be the how around the cylinder at Re = 2. 
They considered Reynolds numbers of 100 and 1000 and chose to hnd either the optimal 
constant rotation rate or sinusoidal rotation amplitude of the cylinder, over horizons of up 
to 5 time-units (i.e. roughly one vortex shedding period). The optimal sinusoidal forcing 
parameters (0 = 3.25, St = 1.13, for Re = 100) are diherent from the ones found by He 
et alM but of the same order of magnitude. Protas et alM also studied this problem but 
considered a free waveform for the control and minimized the sum of the control power 
and the power required to overcome drag, for Reynolds numbers of 75 and 150. They 
considered horizons of up to roughly one vortex shedding period and reported not obtaining 
signihcant improvements by further increasing the length of the horizons. The resulting 
control waveform is discontinuous between each horizon so no clear physical mechanism was 
extracted to explain how the rotation affects the flow-held. Nevertheless, the drag is reduced 
by 7% at Re = 75 and 15% at Re = 150, and with a positive energy balance. On the other 
hand, the amplitude of the control waveform does not seem to be changing signihcantly over 
time and the tangential cylinder velocity is typically of the order of 0.05 to 0.1 for Re = 75 
and of 0.1 to 0.2 for Re = 150. Finally, Bergmann et alr^^ used reduced order models 
based on a Proper Orthogonal Decomposition of the how-held to study the same problem. 
They chose to minimize the turbulent kinetic energy in the how-held and the mean drag of 
the cylinder at a Reynolds number of 200. As the cylinder rotation was constrained to be 
sinusoidal, the control converged towards the energy inefficient open-loop optimum. 

In the present investigation, a similar setup as the one in Protas et is selected, but 
signihcantly longer control horizons than in all known earlier studies are considered (up to 
100 convective time-units). In Sec. HIl we introduce the numerical approach and problem 
setup used for the adjoint-based optimizations. In Sec. IIIIl we show the ehectiveness of 
simply increasing the horizon length by choosing the same Reynolds number {Re = 75) 
and cost function (based on total required power) as Protas et al.—, but with horizons of 


4 


50 convective time-units. We then show that convincing results can also be obtained with 
different cost functions by minimizing squared lift or squared drag with an unconstrained 
control waveform, at a range of Reynolds numbers. The physical mechanism behind the 
performance of the optimal control is then analyzed. Finally, the influence of the cost 
function, of the length of the optimization horizon (10 to 100 convective time-units), and of 
the Reynolds number (between 75 and 200) are all discussed. 


II. PROBLEM FORMULATION AND NUMERICAL METHOD 

The two-dimensional incompressible Navier-Stokes simulations presented in this article 
were run using a finite-volume immersed-boundary fractional step algorithm developed by 
Taira and Colonius^>^. For this study, adjoint solvers were developed based on earlier work 
of Joe et al.—’^ and Ahuja & Rowley^. In this section we provide an overview of the setup 
and procedures used and include further details about the solvers in Appendix and about 
the adjoint equations in Appendix |Bl 

A nominal set of parameters was used in all simulations presented in this article. The 
adequacy of these parameters was therefore checked as summarized in Table [H Here, the 
vortex-shedding Strouhal number and mean drag coefficient of a stationary two-dimensional 
circular cylinder in a free-stream with a Reynolds number between 75 and 200 were computed 
with the nominal parameters, as well as with a more refined set of parameters, and the results 
were compared to those found in the literature. A good convergence and agreement with 
previous studies is obtained with the nominal set of parameters for both Reynolds numbers. 
For the nominal (refined) runs the time-step is 0.002 (0.001) and the flow is solved on 8 
nested grids of sequentially increasing size and correspondingly decreasing resolution (as 
introduced in Appendix IA|) . Each grid is composed of 240 x 120 (600 x 300) cells and the 
extent of the smallest grid is 4.8 x 2.4 (6.0 x 3.0) cylinder diameters. A plot of the cylinder’s 
location within the first 3 (out of 8) nested grids levels in shown in Fig. [H As the adjoint 
solution is advected in the opposite direction to the forward one, the cylinder is located 
at the center of each grid level. The immersed-boundary points that define the cylinder is 
shown in Fig. [21 where it is superimposed on part of the mesh of the finest grid level. 

The adjoint-based optimization procedure considered in the present work aims to mini¬ 
mize a cost function of the form J'{x, (j)) = X{x, which depends on the state of the 
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FIG. 1. Location of the cylinder within the first 3 nested grid levels (out of 8). 
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FIG. 2. (Color Online) Uniform Cartesian mesh of the finest grid level, in the vicinity of the 
immersed-boundary points (red dots) used to define the cylinder surface (of which a quarter is 
shown here). 


flow x{t) and control waveform 0(f) as defined in Appendix |Bl In this case, the integral of 
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TABLE I. Strouhal number of vortex shedding and mean drag coefficient of the unforced flow 
around a stationary circular cylinder with the nominal and refined set of parameters, compared 
with values obtained in previous studies. 


Simulation 

Re 

Strouhal Number 

Mean Cd 

Nominal 

100 

0.163 

1.330 

Refined 

100 

0.164 

1.327 

Braza et al.— 

100 

0.16 

1.36 

Henderson^i 

100 

0.166 

1.350 

Park et alM 

100 

0.164 

1.33 

He et alM- 

100 

0.167 

1.353 

Nominal 

200 

0.194 

1.327 

Refined 

200 

0.196 

1.356 

Braza et al^ 

200 

0.20 

1.39 

Henderson^ 

200 

0.197 

1.341 

He et ali^ 

200 

0.198 

1.356 

Bergmann et 

200 

0.200 

1.390 


the power, the squared lift or the squared drag over the control horizon were used: 


Jv — 

{ (Power)+ (Power)Control } 

( 1 ) 

Jc = 

j \ (LJtf dt, 

( 2 ) 

Jv = 

rT ^ 

J 2 dt, 

( 3 ) 


The cost function is minimized using a single degree of freedom but time-dependent control 
4>{t): The tangential velocity of the cylinder surface. Note that there is no control penaliza¬ 
tion term in (|2]) and ([3]). Although in general there is no guarantee that the problem will 
be well-posed without regularizing the control effort, it was found that for the setup pre¬ 
sented here, no such term was necessary for long enough horizons (50 convective time-units 
or more), in order to obtain a converged optimal solution with a bounded control effort. 

In order to locally minimize these cost functions, the gradient of the cost with respect 
to the controls is determined in the standard manner— as detailed in Appendix [Bl Using 
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calculus of variations, this procedure yields a set of conditions that must be satished in order 
to identify the optimal control. First, the spatially discretized version of the vorticity form 
of the Navier-Stokes equations dl]) must hold throughout the simulation: 

{ (h = V X (u X (u) — Re~^'V x (V x cu) + V x /, 

Ub = 4>{t)ut. 

Here u refers to velocity, a; = V x u to vorticity and co = doj/dt, ub is the velocity on the 
body surface, and (j){t) is the control, which as mentioned above imposes a tangential velocity 
along the body surface, mimicking cylinder rotation: ut is a locally dehned unit vector, 
which points in the direction that is tangential to the body surface in the counterclockwise 
direction. The immersed body force / is used to impose the jet velocity along the body 
surface. 

Second, the corresponding adjoint Navier-Stokes equations must also hold throughout 
the simulation. These can be shown to be a spatially discretized version of the following 
equations (further details in Appendix iBl): 


-Cj+ = V X (o; X n+) - (m+ x m) - x (V x a;+) + V x /+, 

1 +_ + 

where refers to adjoint quantities, is the adjoint “slip velocity” at the body surface. 
The adjoint equations are forced through this term, and this forcing depends on the cost 
function (as again shown in Appendix iBll. 

Finally, the gradient of the augmented cost function with respect to the controls Q must 
be 0. The gradient can readily be computed given the data from the forward and adjoint 
simulations. However, (f){t) is in general not optimal so Q is in general non-zero. It is instead 
used to iteratively update the control waveform using a conjugate gradient approach. 

The full adjoint optimization procedure can thus be outlined as follows: 


1. A cost function J'{x, 0), an initial control guess 0(t), and a starting condition a;(0) are 
dehned; 

2. The Navier-Stokes equations are solved in order to evaluate x{t) from t = 0 to t = T, 
starting from the initial condition a:(0), and with the current control guess (p{t). The 
corresponding cost J'{x,(j)) is also evaluated; 
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3. The adjoint equations are solved backwards in time, in order to evaluate x~^{t) from 
t = T to t = 0, with the “initial” condition x~^{T) = 0. The control gradient Q is 
calculated from the results of the forward and adjoint simulations; 

4. The optimal control update distance is calculated iteratively based on the calculated 
gradient Q{x,x ~^using a line minimization algorithm (more details below); 

5. The control guess 0(t) is updated according to the results of the line minimization, 
allowing the new cost to be evaluated; 

6. Steps 2. to 5. are repeated using the updated control guess 0(t) until convergence; 

7. The starting time is advanced by about 2/3 of the control horizon (in this case) and 
the state of the controlled flow-held at that time is used as the new initial condition. 
The last third of the converged optimal control signal is discarded. The new initial 
control guess is set to zero (more details below); 

8. Steps 1. to 7. are repeated for the next horizon. 

More details about this procedure are included for instance 

The line minimization algorithm mentioned in the 4th step above is necessary as the 
control gradient only provides information about the direction in which the control needs to 
be updated in order to optimally reduce the cost function (locally). A further inner iteration 
is then required to hnd the optimal magnitude of the control update in the gradient direction. 
A steepest gradient algorithm can be used, but it is more efficient to use a conjugate gradient 
approach^. Typically, a Brent line minimization is then used to find the optimal update 
distance^. However, in order to make this searching procedure more efficient, a generalized 
version of Brent’s search algorithm was developed and used in the present work, whereby 
several forward simulations are run together in parallel. 

In order to obtain optimization results that are longer than one horizon, a new initial 
control guess and starting condition are necessary to start the optimization at each horizon. 
As mentioned in Step 7. above, the initial condition for each horizon i of total length Tj 
was chosen to be a snapshot of the converged solution of the previous horizon i — 1 such 
that fj(0) ~ (2/3)Tj_i. A large enough fraction of the optimal control solution should be 
discarded to ensure transients from the adjoint simulation do not affect the retained solution. 
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Here, keeping about 2/3 of the optimal control was chosen as an adequate value to truncate 
the part of the solution that was most affected by adjoint transients. Reasonably small 
variations in this value would only be expected to change the performance of the control 
slightly. As it was chosen to set the hrst control guess to zero rotation throughout each 
new horizon, an initial discontinuity is to be expected in the optimal control signal since 
0((2/3)Tj_i) 7^ 0 in general. 

The gradients computed with this code are not the exact discrete gradients of the forward 
problem. Indeed, the time-stepping scheme is not exactly self-adjoint, as the base-flow is 
unsteady and the same multi-step marching method is used for the forward and adjoint 
solvers, as described in Appendix Additionally, the same nested grid procedure was used 
in the forward and adjoint solvers, but this procedure is in fact not self-adjoint either, as 
noted by Ahuja and Rowley^. Finally, solving the adjoint equations requires storing the 
entire solution of the forward problem (the unsteady base-flow) since u and uj appear in the 
advection terms of Eq. ([5]). To reduce the computational requirements of this procedure, 
we only store every 10th snapshot of the forward solution and use linear interpolation to 
reconstruct it at a cheaper cost. The error associated with this procedure is very small, as 
the base-flow evolves at a much slower rate than the time-step. For instance, at Re = 100, 
using interpolation introduces a relative error in the gradient that is typically of the order 
of 0.001% to 0.005% (with a control guess of = 0 and a horizon of 10 convective time- 
units). Nevertheless, as a consequence of the points mentioned above, we can only expect 
to obtain an approximation to the exact discrete optimal solution of the problem. 

In order to check that the gradients obtained with the code are of an acceptable accuracy, 
a standard finite-difference check— was performed on the adjoint code: if the forward and 
adjoint equations are satisfied, then the following Taylor expansion of the cost function can 
be written: 

J{(p) - J{(p + 54>) = [ + 0{6(j)‘^)) dt, (6) 

Jo 

where Q{t) is the control gradient. The nominal control is chosen to be 0(t) = 0 and 
S(j){t) is an arbitrary control perturbation, chosen to be a sinusoidal function of the form 
6(j){t) = esin {{27rSt) t), for small enough e. Here e = 10“® was used and Strouhal numbers in 
the range 0.1 < St < 0.5 were tested. The error between the two sides of Eq. ([6]) normalized 
by \J{.P) — + d(j))\ was checked and found to be typically of the order of 0.5%, 1% and 
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TABLE II. Summary of results from the optimizations: the horizon length is given in convective 
time-units, J refers to the cost functions defined in Eq. (EEIED. Re is the Reynolds number, Co 
and Cl are the drag and lift coefficients, and cf) is the control amplitude. The values quoted here 
refer to the reduction in mean Cd and the relative change va. Cl amplitude (with respect to the 
unforced case at the same Re), as well as the amplitude of the tangential velocity cj). All values 
are estimated towards the end of the last considered horizon for each Run, ignoring the part of 
the response where adjoint transients become significant. The results are also compared to those 
obtained by Protas et al^. 


Run 

Horizon 

J 

Re 

Cd 

Cl 

<P 

Protas et alM 

6 

Jr 

75 

7% 

- 

0.05 - 0.1 

Protas et al.^ 

6 

Jr 

150 

15% 

- 

0.1 -0.2 

1 

2x50 

Jr 

75 

13% 

-99% 

0.0007 

2 

2x50 

Jv 

75 

13% 

-95% 

0.002 

3 

6x10 

Jv 

100 

> 60% 

> +1900% 

> 2.0 

4 

50 

Jv 

100 

16% 

-89% 

0.04 

5 

2x50 

Jv 

100 

19% 

-98% 

0.001 

6 

100 

Jv 

100 

10% 

-69% 

0.10 

7 

50 

Jc 

100 

6% 

-73% 

0.07 

8 

50 

Jv 

200 

9% 

-58% 

0.21 

9 

2x50 

Jv 

200 

9% 

-46% 

0.14 

10 

100 

Jv 

200 

7% 

-38% 

0.12 

11 

50 

Jc 

200 

11% 

-64% 

0.21 


3% or less for T = 10, T = 50, and T = 100, respectively for all three cost functions given 
in Eq. (IH Ej |3]) and at Re = 100. Due to the fact that the error seems to increase with T, 
we can expect to obtain more suboptimal results with longer horizons. 

III. RESULTS AND DISCUSSION 

Several optimizations were run in order to investigate the influence of changing the cost 
function, the horizon length, and the Reynolds number, as summarized in Table [TTl The 
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FIG. 3. Results from Run 5. (a) Unforced wake vorticity contours, (b) Optimally controlled wake 
vorticity contours. Solid (dashed) lines represent positive (negative) vorticity. Contours shown for 
vorticity values of -2 to -0.5 and 0.5 to 2 in increments of 0.5. The vortex shedding is almost fully 
suppressed in the controlled case. 

performance of the optimizations is qnantified by the amonnt drag redaction and snppression 
of lift flnctnations obtained towards the end of the optimized simnlations. 

As shown in Rnns 1 and 2, in which the Reynolds nnmber was chosen to be the same as 
the one nsed by Protas et al.— {Re = 75), improved resnlts are obtained just by increasing 
the length of the control horizon, both when the same cost function is used (Run 1 with 
J'p) and when the drag cost function is used (Run 2 with jTp). Table UTI shows that only a 
small amount of control is required towards the end of the horizon, and the lift fluctuations 
and mean drag are both reduced by a larger amount than with the short horizons used in 
previous work. 

In Sec. IIII Al we investigate how the optimal control leads to such large reductions in 
the mean drag and lift fluctuations, by focusing on Run 5, where the Reynolds number is 
increased to Re = 100 and JT” = We then discuss the influence of the cost function, 
the Reynolds number and the horizon length on the overall optimization performance in 

Sec. Jirm 

A. Suppression of vortex shedding 

In Fig. E] the unforced flow-held (a) is compared to the optimally controlled how-held 
(b) of Run 5. The intensity of the vortex shedding in Fig. [3](b) is reduced and the how 
appears to be much more symmetric. The optimal control signal is shown in Fig. 111(a), 
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while the optimally controlled drag and lift coefficient signals are shown in Fig. IH^b) and 
(c) respectively. In Fig. Ht^a), the drag of the unstable (steady) equilibrium is also plotted 
(this value was computed by using Selective Frequency Damping^ii^). Clearly, the controlled 
flow-field has been stabilized to a drag state that is close to that of this unstable equilibrium. 

From Fig. [3l Fig. 01 and the values from Table m it is clear that in Run 5 the vortex 
shedding is almost fully suppressed and that the lift oscillations and the mean drag are 
both reduced significantly. Moreover, towards the end of the second optimization window, 
the tangential velocity on the cylinder surface is of the order of 0.1% of the incoming flow 
velocity. Comparing these results to previous work with similar setups, only Protas et alM 
obtained comparable drag reductions without needing an excessively large amount cylinder 
rotation. However, the control waveform obtained by Protas et ai— is discontinuous due 
to the chosen setup for the optimizations: only short optimization horizons were used. In 
the present case, almost no cylinder rotation is required to keep the Von Karman wake 
suppressed for large times. Moreover, we obtain a smooth control waveform, which clearly 
has a main frequency component that corresponds to the vortex shedding frequency. We 
therefore proceed to investigate in what way the wake is affected by the control. 

In previous optimal control studies^^'^S.'^^, the aggressive optimal control obtained can 
be seen as essentially an open-loop strategy with optimized (constant) parameters: purely 
harmonic forcing cannot suppress the wake instability, it can only alter it to a more favorable 
drag state with suitably large rotational velocities. In order to test whether the low drag state 
reached with optimal control can also be obtained with similar open-loop control signals, two 
open-loop control waveforms that approximate the optimal control were designed, as shown 
in Fig. m the first is a sinusoidal signal whose frequency, amplitude, and phase approximately 
match the optimal control at f = 0 and the second is an exponentially decaying sinusoidal 
waveform (also with approximately matching frequency, amplitude and phase at t = 0). 
Upon applying these control signals to the cylinder, the drag is initially reduced in both cases, 
but it returns to its unforced value for the exponentially decaying sinusoidal case. In the 
sinusoidal case the control actually increases the drag by 9% for large times, with respect to 
the unforced case. This type of behavior is common: the vortex shedding frequency naturally 
locks-in to the forcing frequency when it is close to the unforced shedding frequency, thereby 
amplifying the vortex shedding intensity (e.g.— ). 

These simple tests therefore seem to confirm that there is an intrinsic need for feedback 
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FIG. 4. Results from Run 5. (a) Optimal control waveform and resulting drag (b) and lift (c) 
coefficients, showing a clear suppression of the lift oscillations, a significant reduction in the mean 
drag, and a control waveform that tends to nearly zero amplitude for large times. In (b) the drag 
of the unstable equilibrium of the flow-field, computed by using Selective Frequency Damping2ii2^, 
is also plotted as a reference (dashed line). 

information to be available to the control for it to damp the vortex shedding effectively. 
Thiriaii^ studied the influence of the rotation of the cylinder on the wake vorticity in the 
open-loop case and argued that when the optimal drag reduction forcing parameters are 
chosen, the rotation of the cylinder has a “destructive” influence on the formation of vortices, 
since the cylinder rotation creates vorticity of the opposite sign to that of the vortex being 
formed. At the same time, the rotation promotes the creation of a vortex in the opposite 
shear layer, thus effectively enhancing the synchronization of the shedding of vortices in 
the wake. On the other hand, when the rotation is locked-in with the vortex shedding, 
the interaction can be seen as “constructive” since the cylinder rotation creates vorticity of 
the same sign as that of the vortex being formed and reduces the amount of vorticity in 
the opposite shear layer, thus aggravating the shear layer interaction and strength of the 
vortices. 
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FIG. 5. (Color Online) Results using harmonic approximations to the optimal control waveform of 
Run 5 (also run with Re = 100): (a) Control signals and resulting (b) drag and (c) lift coefficients. 
Black solid line: optimal control. Red dash-dotted line: exponentially damped sinusoidal forcing. 
Blue dashed line: sinusoidal forcing. Despite a similar initial behavior, neither of the open-loop 
approximations are able to suppress the shedding sustainably. 

Comparing the direction of the rotation with the vortex shedding phase for the open- 
loop sinusoidal case described above (blue dashed line in Fig. |5]) the same conclusions can be 
drawn: as shown in Fig.[6](a) and (c), for t < 10, the rotation is in the “destructive” direction, 
while for t > 60, (Fig. [6](b) and (d)) the vortex shedding is locked-in and the rotation is in 
the “constructive” direction. This suggests that matching the cylinder’s rotation direction 
to the phase of the vortex shedding is a promising strategy to suppress the wake fluctuations 
without requiring the full optimal control framework. 

One way to interpret this conclusion is that the mechanism leading to the suppression 
of vortex shedding in the open-loop experiments of Tokumaru and Dimotakis^ and the 
subsequent related studies is effectively the same as the one allowing suppression of vortex 
shedding in the energy-efficient optimal control case shown here. In the open-loop case, a 
large amount of input energy is required in order for the interaction between the destructive 
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FIG. 6. Vorticity contours for the sinusoidal forcing case with amplitude, frequency and phase 
matching the optimal control for early times. Solid (dashed) lines represent positive (negative) 
vorticity. Contours shown for vorticity values of -2 to -0.5 and 0.5 to 2 in increments of 0.5. (a,b) 
Maximal value of positive (anti-clockwise) rotation. (c,d) Maximal value of negative (clockwise) 
rotation. (a,c) In the first few shedding periods, the cylinder rotation is “destructive” and reduces 
the intensity of the shedding. (b,d) Once lock-in occurs, the rotation becomes “constructive” and 
reinforces the shedding. 

sinusoidal rotation and the vortex shedding to become stable. On the other hand, similar 
performance is possible with a much smaller control effort if information is fed back to 
the controller as this allows the unstable destructive interaction to be maintained without 
requiring the fast high-amplitude rotation that is necessary with open-loop forcing. 

B. Impact of optimization and simulation parameters 

The results in Table m show that at a Reynolds number of 100, minimizing the RMS 
drag is more efficient in this case than minimizing the RMS lift (see Run 4 compared to 
Run 7). For a Reynolds number of 200 however, the performance is comparable. For both 
cost functions and Reynolds numbers, the cylinder rotation is slow, of small amplitude 
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and roughly at the vortex shedding frequency, and the drag is signihcantly reduced due to 
the suppressed intensity of vortices in the wake. In previous adjoint optimization studies 
(e.g.— d2.) it was also found that it is not always immediately obvious why one cost function 
performs better than another one, due to the highly non-convex and nonlinear nature of the 
problem. Clearly here, we expect both Jn and jTb to be minimized if the vortex shedding is 
somehow completely suppressed, but there is no guarantee that using the two cost functions 
will in fact lead to similar modihcations in the flow-held, as the path to the local minimum 
reached and even the minimum itself may be diherent. The present results therefore bring a 
further conhrmation that it is worth running optimizations with several cost functions even 
if all of them are expected to be reduced by the desired how behavior. Note that although 
all simulations were run for enough iterations to appear converged, it is possible that further 
reductions in the cost could be obtained by computing a larger number of iterations, since 
changing J also ahects the convergence properties of the optimization problem. 

Table [ITl also shows that increasing the Reynolds number from Re = 100 to Re = 200 
is overall detrimental to the relative performance of the control. From linear stability 
analysis^, the real part of the unstable eigenvalues associated with the vortex shedding 
grows with the Reynolds number and several studiesl^i^ also reported a reduction of closed- 
loop performance as the Reynolds number is increased. In the present case, the two shear 
layers are more unstable for Re = 200 than Re = 100, so even if the cylinder rotation is 
able to initially enhance the symmetry of the wake and delay the interaction between the 
shear layers, the instability is not as readily suppressed for the entire wake region and the 
drag reduction is therefore not as signihcant. At both Re = 75 and Re = 100, the flow can 
be considered to be effectively stabilized (as suggested by Table [TTl which shows that the 
lift fluctuations are almost fully suppressed). At Re = 75 however, as the Reynolds number 
is lower than at Re = 100, not only is the vortex shedding is less intense, but the unstable 
equilibrium also has a higher Crr^. As a result, a smaller relative drag reduction is possible 
by stabilizing the wake in this case. 

A key difference between previous work and the setup chosen here is that in this study, 
it was chosen to use much longer optimization horizons. Protas et al— found that the 
horizon length should be longer than one vortex shedding period but only obtained marginal 
improvements by further increasing it. Bewley et alM argued that if the horizon is too long 
the optimization can become excessively non-convex and computationally expensive, leading 
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FIG. 7. Results from Run 3. (a) Optimal control waveform and resulting drag (b) and lift (c) 
coefficients. The optimal control leads to a non-zero-mean, fast rotation optimum. 

to suboptimal performance in practice. Furthermore, inaccuracies in the computed gradient 
have the potential to accumulate and stall the optimization for very long horizons, leading 
to more suboptimal results. In this article, optimizations with horizon lengths of 10, 50 and 
100 convective time-units were run. 

Considering the 100 convective time-unit optimizations, a similar behavior to the 50 
convective time-unit case can be seen in Table [TTl but with an inferior performance (see Run 
5 compared to Run 6 and Run 9 compared Run 10). This correlates well with the fact that 
100 convective time-units is indeed an excessively long horizon in this case. As mentioned 
above, it is possible that additional iterations may still result in a further reduction in the 
cost. Nevertheless, the present results conhrm that such long horizons can indeed result in 
prohibitively slow convergence of the optimization problem or suboptimal converged results. 

When on the other hand the horizon length is 1 convective time-unit (not shown here), 
i.e. shorter than a vortex shedding period, it was found that the problem does not converge 
at all, as it is ill-posed from the hrst optimization horizon, without a control penalty term. 
With a horizon length of 10 convective time-units however (Run 3), the solution of the 
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FIG. 8. Vorticity contours corresponding to Run 3. Solid (dashed) lines represent positive (neg¬ 
ative) vorticity. Contours shown for vorticity values of -2 to -0.5 and 0.5 to 2 in increments of 
0.5. The flow-field is strongly asymmetric due to the large nearly constant rotation of the cylinder, 
which leads to reduced drag. 

optimization tends towards a well-known low-drag flow-state with nearly constant rotation, 
as shown in Fig. [71 The control waveform has discontinuities at the start of each horizon 
but the general behaviour of the control is still clearly visible. Homescu et alM and Kang 
et ai— also showed that the cylinder drag can be reduced signihcantly and the vortex 
shedding suppressed completely by using sufficiently fast constant rotation of the cylinder 
in one direction. Unsurprisingly, this results in a signihcant constant lift force on one side of 
the cylinder, as shown in Fig. [71(c). Both studies found an optimal rotation rate at Re = 100 
of 0 ~ 1.85. Here the rotation rate quickly increases until it reaches a value that is close 
to this “optimum”, but then continues to slowly increase. After 40 convective time-units, 
the rotation rate is 0 ~ 2.1 and the cylinder is still slowly accelerating. It is expected 
that computing further horizons would lead to a continued increase in this quasi-steady 
rotation rate, since the drag has been shown to decrease further for even higher rotation 
rates^. The physical mechanism reducing drag and stabilizing the wake in Run 3 is therefore 
effectively an open-loop strategy, which does not require closed-loop control but does require 
a signihcant amount of input energy. 

These results show that the choice of horizon length can have a crucial impact on the 
solution of the optimization and suggest the existence of another important timescale than 
the vortex shedding frequency. It seems reasonable to expect a control that efficiently 
suppresses vortex shedding to do so over a length of time that is longer than one vortex 
shedding period. For instance a control waveform that eventually fully stabilizes the wake 
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will not necessarily lead to the maximum drag reduction in the hrst few convective time- 
units. 

One way to visualize this time scale is to realize that the adjoint flow-held will be a tran¬ 
sient simulation, even if based on a fully periodic limit-cycling base-how. From a numerical 
point of view, given that the gradient of the controls is in general dependent on both the 
forward and adjoint state, the horizon should be long enough for its value at the start of 
the simulation not to be signihcantly ahected by any further increase in the horizon length. 
Comparing the forces of the adjoint simulation from the hrst iteration of Run 4 to those 
of Run 6 (Fig. [9]) it is clear that the adjoint simulation experiences transients for longer 
than one vortex-shedding period before the two sets of forces overlap (recall that the adjoint 
simulation runs backwards in time). 

For the last few vortex-shedding periods of the forward simulation, the control will be 
therefore optimized according to what can be considered a “short term” control strategy, 
which may not be sustainable. In some cases, having several short horizons might just lead 
to a suboptimal but similar control waveform to the long horizon case, but in others, the 
two physical mechanisms leading to a reduction in the cost might be different altogether 
(as with the 10 convective time-unit horizons in Run 3). Note however that as the control 
can simply tend towards two different local minima, it is not a priori obvious which one will 
have the lowest cost (indeed, the drag is reduced more with fast constant rotation than with 
the oscillatory rotation in this case). 


IV. CONCLUSIONS 

In this article, optimal control is applied to the flow over a cylinder, where the total 
power, the mean-squared lift or the mean-squared drag were minimized using unconstrained 
rotation of the cylinder. The optimization horizons were chosen to be longer than in previous 
studies (where they typically are at most of the order of a vortex shedding period). 

It was shown that at Reynolds numbers between 75 and 200, cylinder rotation can ef¬ 
fectively suppress vortex shedding in the wake. The optimal control waveform was found 
to be phase-locked with vortex shedding, and the rotation is applied in a direction that 
weakens the vortices shed in the wake. The stabilizing action of the control on the wake is 
therefore similar to the high-amplitude and high-frequency open-loop optimum that is well 
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FIG. 9. Time history of the adjoint force coefficients for the first optimization iteration (no cylinder 
rotation in the corresponding forward run, i.e. 4>{t) = 0) of Run 4 (solid line) and Run 6 (dashed 
line), (a) Adjoint drag C^, (b) Adjoint lift C'^. Note that the adjoint force coefficients are defined 
in an analogous manner to the standard force coefficients (see Appendix!^, but based on the 
adjoint forces at the immersed-boundary points instead. The transients of the adjoint simulations 
can be seen to correspond to several vortex-shedding periods before the two curves start to overlap. 

documented in the literature, but feedback vastly reduces the required rotation rate. 

It was found that this phase-locking behavior can only be achieved in a closed-loop control 
setting. Applying an almost identical but harmonic (i.e. open-loop) forcing signal instead 
results in a drift between the forcing and shedding phase, and eventually increases drag. 
On the other hand, the optimal forcing decreases the drag by up to 19% compared to the 
unforced case and the amplitude of the actuation required to keep the flow in this stabilized 
state approaches zero, with tangential velocities on the cylinder surface of the order of 0.1% 
of the inflow velocity towards the end of the considered optimization window. This shows the 
importance of keeping the control waveform fully unconstrained, as previous studies where 
only periodic solutions were considered reached the high-amplitude and frequency open-loop 
optimum rather than the stabilizing control obtained here. 

The analysis of the influence of different optimization parameters conhrmed the fact that 
the setup of the optimization - including the choice of cost function and horizon length - 
can strongly affect the final results. In particular, changing the horizon length can even 
change the local minimum that is identihed and hence the mechanism through which the 
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cost function is minimized. In this case, it was found an optimization window longer than 
the vortex shedding frequency was required in order for the converged optimal control to 
stabilize the wake in an efficient manner. It was found that the transients of the adjoint 
simulation can provide some important information regarding the adequacy of the chosen 
horizon length. If the horizon is long enough, the state of the adjoint solution and hence the 
control gradients will not be strongly affected by any further increase of the horizon length. 
On the other hand, if the horizon is too short, this can result in a “short term” optimization, 
whereby the reduction in the cost must be measured quickly but may be associated with a 
different control mechanism to the one corresponding to the “long term” optimum. In this 
case, with horizons of 10 convective time-units, the drag was minimized with a fast, nearly 
constant rotation of the cylinder, whereas with longer horizons of 50 convective time-units 
or more, the low amplitude, approximately zero-mean, oscillatory rotation of the cylinder 
described above was obtained. 

In this article, we have demonstrated the importance of choosing a long enough optimiza¬ 
tion horizon when considering adjoint-based optimal control problems, especially when they 
are based on unsteady nonlinear base-flows. We have shown that increasing the horizon 
length not only leads to smoother waveforms but also to potentially large improvements in 
the results. Furthermore the smooth control waveforms make it possible to interpret the 
physical mechanisms that are responsible for the effectiveness of the control. Our work also 
confirmed that excessively long horizons can be detrimental to the convergence of the re¬ 
sults and that changing the horizon length can result in a control waveform that uses an 
altogether different mechanism thus reaching a separate local minimum in the cost function. 


Appendix A: Forward and adjoint solvers and procedures 


The immersed-boundary fractional step solver developed by Taira and Colonius^Ii^S solves 
the vorticity form of the incompressible Navier-Stokes equations (jl]) numerically. Further 
information regarding the solvers can also be found for instance in4 and^. 

Using a uniform staggered Cartesian mesh in regions close to the cylinder, a Dirichlet 
boundary condition is imposed at a set of arbitrarily defined Lagrangian points, which define 
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(Al) 


the body surface, as shown in Fig. [2l The spatially discretized form of the equation reads: 

7 + C^E^f = + C^Af{q, 7 ) + be, 

ECs = A(f)Ut. 

Here, 7 represents the discretized circulation around any given cell (evaluated at the cell 
center) and s is the discretized streamfunction (also evaluated at the cell center), ut is a 
unit vector dehned at each immersed-boundary point, which is oriented in the tangential 
direction to the body surface in the counterclockwise direction and 0 is the time-varying 


control amplitude. C is the discrete curl operator, constructed so that q = Cs, where 

1 T 

^ ^ L ^ " J 

at the left and bottom cell edges respectively. Similarly, we have 7 = C^q, as well as 

T 

are the 


Qx Qy whose X and y components q^ and qy are evaluated 


s = {C^ C) ^ 7 . The discrete Laplace operator is given by —C^ C^. f = 


r f 

J X J 1 


y 


discretized immersed forces with a scaling factor, evaluated at the immersed body points and 
fx and fy are the x and y components of the vector at each point. E and E"’" are interpolation 
and regularization operators that allow evaluating the quantities dehned on the Cartesian 
mesh at the immersed body point locations and vice-versa. Af{q, 7 ) is the discretization of 
the nonlinear advection operator u x uj. Finally, be are the boundary conditions imposed at 
the edges of each domain and /3 = l/(A^i?e), where A is the uniform grid spacing and Re 
is the Reynolds number. 

In the present work, the body is moving horizontally with a velocity Uoo through the huid. 
Staying in the body reference frame, the (potential) hux of the incoming how is therefore 
—goo and the total hux vector thus becomes q — q^o- Equations flAip therefore become: 

7 + C^E^f = -(3C^C^ + C^M{q - goo, 7 ) + be, 

E{Cs - goo) = A(t)Ut, 
which we can write: 


(A 2 ) 


(A3) 


7 = 

E{Cs - goo) = A(j)Ut. 

A second-order (explicit) Adam-Bashforth scheme is used to discretize the nonlinear term 
and the (implicit) Crank-Nicolson method is chosen for the linear terms. The time-stepping 
scheme is based on a diagonalization of the Laplacian using the fast Sine transform and a 
fractional step method, where the Poisson equation to be solved is only of the dimension of 
the number of immersed body force vector /. 
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In order to use far-field boundary conditions, an excessively large uniform Cartesian 
domain would be required for most open flows. Instead, the flow is solved on a series of 
nested uniform Cartesian grids, where each grid level is identical to the previous one, but 
has twice the physical extent and hence half the resolution, as illustrated in Fig. [TJ The 
coarsest grid is chosen to be large enough for the far-held boundary conditions to be justihed 
there. Boundary conditions for all remaining grid levels can then be interpolated from the 
next grid level. The immersed-boundary points are only assumed to exist on the hnest grid. 
At each time-step, the circulation is thus hrst obtained there, allowing the solution to be 
coarsihed onto all other grid levels. 

It is straightforward to evaluate the total force acting on the body in a given direction 
(and therefore the corresponding force coefficient) for a given how state. This force can be 
obtained by summing the components of the forces at all the immersed-boundary points in 
the direction of interest. As a result, the three cost functions we consider in this article, 
introduced in Eq. ([H El |3]), are calculated in the following way: 

Jv = 

Jv = ^fx^fxdt. 

As shown in Appendix |Bl the discretized adjoint equations (1B3|1 are similar in form to 
the forward equations flAlll . except for the advection terms and the fact that the adjoint 
equations run backwards in time from t = T to t = 0. A very similar solver to the one 
introduced here is therefore used for the adjoint equations, including the discretization and 
time-stepping schemes, the far-held boundary conditions, and the nested-grid algorithm. 


Appendix B: Adjoint Eqnations 

Starting from the spatially discrete but time-continuous forward equations flA3l) . a stan¬ 
dard adjoint-based optimization procedure is used in this article. We hrst dehne a cost 
function that depends on the state variables x and the imposed control (j). In our case, it 
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takes the form: 


J{^A) 


T{x, (f))dt. 


(Bl) 


The state variable comprises of the circulation and the immersed-boundary forces x = 


1 T 


The cost is augmented by Lagrange multipliers (the adjoint state) x~^ = 
constrain the state variables to respect the forward equations 

i-T 


7 


+T j+T 


7 ^ r 

r 

that 


J {x, 0, x^) = / |x(a;, 0) - 7 +"^ (C^C) ^ [7 - J^{x)\ - /+'^ [E{Cs - goo) - A0hi] jdf, 


'H{x,(j),x~^)dt. (B2) 

Note that for the momentum equation, we use a non-trivial symmetric inner-product matrix 
so that ( 7 , 7 ) = 7 ’^(C^C )“^7 = q^q, which is directly related to the kinetic energy 
of the flow. All variations of Eq. flB2p with respect to x, x^ and 0 must be zero for the 
solution to be optimal. Thus, setting variations with respect to x'^ to zero hrst, we obtain: 



I i = 

y E{Cs - goo) = 

In other words, the forward equations must be enforced throughout the horizon. Next, 
setting variations with respect to x to zero can be shown to result in the following adjoint 
equations: 

' - 7 + + = -/?C^C 7 + + (C^C) q* + (C^C) + bc+, 

< ECs+^i^af, (B3) 

_ 7+(r) = (C^C) {^Y, 

where q'^ = Cs~^ are the adjoint fluxes, s’*" = is the adjoint streamfunction, and 

fee"*" are the adjoint boundary conditions. All quantities here are dehned in an analogous 
manner to their counterpart in the forward equations. Equations flB3D must also be satished 
throughout the horizon, but they are integrated backwards in time from the initial condition 
at t = T (given by the third equation in flBSp h For the cost functions considered here, we 
have simply 7 ’'"(T) = 0, since we do not have a cost on the hnal flow condition. In fact, 
it can be shown that these equations can also be obtained by directly discretizing Eq. (|5]), 
which are their continuous equivalent. 
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Note that the adjoint equations are very similar to the forward equations flAlj) . except 
for the negative time derivative and the advection term. This term can be shown to be 
a discretization of V x (cj x m+) — (m"*" x u) and therefore depends on the (unsteady, 

nonlinear) state of the forward simulation at time t (since u and u are needed to compute 
it). The entire forward base-flow 7(f) for all t such that 0 < t < T is thus required to 
compute the adjoint solution. 

Note also that the second equation in fIBSp sets the “slip flux”: ECs^ = = (^dX/df'^ 

at the immersed-boundary points and is dependent on the cost function. For the three cost 
functions considered in this article - i.e. J7^p, J7/;, and JT©, given in Eq. firiEll^ - the slip flux 
is given respectively by: 

= ^(put + AEqoo, 

1 T 


^slip,C 

^slip^Ti 


0 

/J 0 


The right-hand-side forcing term of the adjoint momentum equation (the hrst equation in 
flBSp i is also dependent on the cost in general: for the cost functions considered in this 
article we have simply dXjd^ = 0. 

Finally, taking variations with respect to 0 yields the control gradient: 




d(\)) 


The gradient will only be zero if the control is exactly optimal, which is not the case in 
general. Its direction is therefore used to update the current control guess in a locally 
optimal manner. For the cost functions considered here - again Jp, JT^, and jTp, given in 
Eq. ffT]|2]l5]l - the expression of the three gradients are respectively: 

Qv = Air 


f+r' 


Qc. = Af-ut, 

Gv = Arr. 
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